%--------------------------------------------------------------------------
%
% computes the leading order (in delta) solution to the basic state and
% then compares to the numerical one. doesn't compare the actual 
% concentrations, however. Instead, compares the average concentration in
% the layer
%
%--------------------------------------------------------------------------

function base_compare_average(p, trans_c)

if nargin == 0
    p = params;
    close all;
end
if nargin < 2
    trans_c = 0;
end
s = base(p);

v = s.y(1:end-1,:);
h = s.y(end,:);
t = s.x;

% change my v into the paper's c
c = v / p.delta / p.beta / (1 - p.beta);

c_mean = (c(end,:) - c(1,:)) ./ h;

h_comp = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);


t_i = logspace(log10(t(2)), log10(t(end)), 40);
h_i = interp1(t, h_comp, t_i);


v_0 = (1 - p.beta) * (1 - 1 ./ h);
va_mean = p.delta * 1/2 * (1 - p.beta) * (-p.beta - v_0) ./ h;
ca_mean = va_mean / p.delta / p.beta / (1 - p.beta);


v1 = sep_var(0);
c1 = v1 / p.delta / p.beta / (1 - p.beta);
c1_mean = (c1(end,:) - c1(1,:)) ./ h_comp;

semilogx(t_i, interp1(t, ca_mean, t_i), 'k*',t_i, interp1(t, c1_mean, t_i), 'ko');
hold on;
semilogx(t, c_mean, 'k', 'linewidth',1);
% semilogx(t_i, interp1(t, ca_mean, t_i), 'k*');
% ylim([1e-5, 1]);
xlabel('$t$', 'interpreter','latex','fontsize',12);
ylabel('$[\bar{c}(\bar{h}(t),t) - \bar{c}(0,t)]/\bar{h}(t)$', 'interpreter','latex','fontsize',12);
if ~trans_c
    l = legend('numerical','$t = O(\delta^{-1})$ asymptotic', '$t = O(1)$ asymptotic', ...
        'location', 'southwest');
else
    l = legend('$t = O(\delta^{-1})$ asymptotic', '$t = O(1)$ asymptotic', ...
        'location', 'southwest');
end
set(l, 'interpreter','latex','fontsize',12);
ylim([-1 0.1]);